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Abstract 

Variational principle for Kolmogorov-Petrovsky-Piskunov (KPP) 
minimal front speeds provides an efficient tool for statistical speed 
analysis, as well as a fast and accurate method for speed compu- 
tation. A variational principle based analysis is carried out on the 
ensemble of KPP speeds through spatially stationary random shear 
flows inside infinite channel domains. In the regime of small root 
mean square (rms) shear amplitude, the enhancement of the ensem- 
ble averaged KPP front speeds is proved to obey the quadratic law 
under certain shear moment conditions. Similarly, in the large rms 
amplitude regime, the enhancement follows the linear law. In particu- 
lar, both laws hold for the Ornstein-Uhlenbeck process in case of two 
dimensional channels. An asymptotic ensemble averaged speed for- 
mula is derived in the small rms regime and is explicit in case of the 
Ornstein-Uhlenbeck process of the shear. Variational principle based 
computation agrees with these analytical findings, and allows further 
study on the speed enhancement distributions as well as the depen- 
dence of enhancement on the shear covariance. Direct simulations in 
the small rms regime suggest quadratic speed enhancement law for 
non-KPP nonlinear ities. 
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1 Introduction 



Front propagation in heterogeneous fluid flows has been an active research 
topic for decades (see [7j , [T7j , [H] , [2H] , [21] , jlH] , [21] and references therein) . 
A fascinating phenomenon is that the large time (large scale) front speed can 
be enhanced due to the presence of multiple scales in fluid flows. Speed char- 
acterizations and enhancement laws have been studied mathematically for 
various flow patterns by analysis of the proto-type models, e.g. the reaction- 
diffusion-advection equations (see 01311111211131111111111312111221121112111201 
and references therein) . The enhancement obeys a quadratic law in the small 
amplitude flow regime, known as the Clavin- Williams relation jjj, which is 
proved to be true for deterministic shear flows [22], [12]; [20]; [21] etc. The 
enhancement is proved to grow linearly in large amplitudes of deterministic 
flows (shear and percolating flows), see [H [HI HU EES 1201 121] etc. However, en- 
hancement exponent 4/3 was reported based on numerical simulation of ran- 
dom Hamilton- Jacobi models (so called G-equation or KPZ model) on fronts 
in weak randomly stirred array of vortices ^S]- Likewise, formal renormal- 
ization group method [2Z] suggested that front speeds may grow sublinearly 
(slower than linear by a logarithmic factor) in the strong random flow regime. 
These two findings raised the issue as to what extent the speed enhancement 
laws in the deterministic flows are valid for random flows. Recently, one 
of the present authors [2H] showed that both the quadratic and linear laws 
hold almost surely for Kolmogorov-Petrovsky-Piskunov (KPP) minimal front 
speeds through white in time spatially Gaussian random shears on the plane. 
Yet the KPP front speeds diverge almost surely in spatially Gaussian ran- 
dom shears. A powerful tool is the variational principle of KPP minimal 
front speeds. 

In this paper, we consider KPP front speeds through random shear flows 
in channel domains D = R x Q, where Q C i? n_1 , n > 2, a bounded simply 
connected domain with smooth boundary. We shall address the enhancement 
laws of the ensemble averaged front speeds. The model equation is: 

u t = A x>y u + B ■ V XjV u+ f(u), (1.1) 

where t G R + , A xy the n-dimensional Laplacian, (x,y) G D. The nonlin- 
earity / = u{l — u), the so called KPP reaction. Other nonlinearities [232| 
will be discussed later. The vector field B = (b(y,u),0) where b(y,uj) is a 
stationary continuous scalar random process in y, its ensemble mean equal 
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to zero. The zero Neumann boundary condition is imposed at dQ: |^ = 0, 
v the unit outward normal. 

For nonnegative initial data approaching zero and one at x infinities 
rapidly enough, the KPP solutions propagate as fronts with speed c* given 
by the variational principle [HJ 123 IS] : 



where /i(A, u>) is the principal eigenvalue with corresponding eigenfunction 
(f) > of the problem: 



L x <f> = A y 0+[A 2 + A%,cu) + f(O)]0 = //(A,u;)0, yeQ, (1.3) 



The variational speed formula (jl.2J) makes possible an analysis of ensem- 
ble averaged random front speeds. Using variational formulas on the principal 
eigenvalue /i, we are able to obtain tight upper and lower bounds on c* and 
estimate E[c*] in terms of moments of suitable norms of the shear over Q. 
For small root mean square (rms) shear, the quadratic enhancement law is 
proved and is explicit in case of the Ornstein-Uhlenbeck (O-U) process when 
n = 2. The linear growth law holds in the large rms regime under weaker 
moment conditions on the shear. In both regimes, the moment conditions 
are satisfied by the O-U process when n = 2. 

The variational formula (J 1.2)) also offers an efficient and accurate way 
to compute a large ensemble of random front speeds. Directly solving the 
original time dependent equation (jl.l|) to reach steady propagating states 
can be both slow and less accurate. Numerical difficulties abound for direct 
simulations due to three large parameters. A large ensemble of random fronts 
requires a large enough truncated channel domain to contain the front over 
large times. Due to occasional random excursions in b, the domain size in x 
has also to be made adaptively large. This can be prohibitively expensive in 
the regime of large rms shears. 

The variational formula (jl.2|) allows us to compute quickly and accurately 
the ensemble averaged speeds in both small and large shear rms regimes 
when n = 2. An interesting difference from the deterministic case is that 
the integral average of b(y,u) in y G Q, i.e. b = b{ui) = J n b(y,u) dy, 
is a random constant not equal to zero. This quantity can be of either sign, 
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and influence greatly the numerical approximation of E[c*] in the small rms 
regime, even though it does not contribute to the exact £?[c*], since E[(b)] = 
0. To assess the speed enhancement accurately, we subtract this random 
constant from each c*(a>) before evaluating the expectation numerically. This 
way, we are able to minimize the errors in approximating E[c*] in a finite 
ensemble. In our computation, b is a discrete O-U process. 

In complete agreement with analysis, we find numerically that the en- 
semble averaged speeds obey a quadratic law in the small rms regime and 
a linear law in the large rms regime. Without the b subtraction technique, 
the computed average speed enhancement in the small rms regime can give 
inaccurate scaling exponents significantly below two. The same technique 
and direct simulations for other nonlinearities (combustion, bistable) suggest 
quadratic speed enhancement in the small rms regime. The computed speed 
ensemble then permits us to study further the enhancement distribution and 
its dependence on variation of shear covariance functions. 

The paper is organized as follows. In section 2, we prove the enhancement 
laws of ensemble averaged speeds based on variational principles, and derive 
a closed form speed asymptotic formula in case of O-U process. In section 
3, we describe numerical methods for computing speed ensemble with the 
variational formula (jl.2J) and speed statistics, then show numerical results. 
We also make comparisons with the prediction of the asymptotic formula and 
with direct simulations. In section 4, we conclude with a remark on future 
works. 

2 Average Speed Asymptotics 

Consider scaling shear amplitude b(y) i— > Sb(y), and denote by c*(S) the 
minimal KPP speed corresponding to shear 5b. Let Co = c*(0) = 2a/ /'(0) 
denote the minimal speed in case of zero advection. If the shear b = b(y) 
has zero integral average over Q, (b) = J n b(y) dy = 0, the corresponding 
minimal speed c*(u) is always enhanced by the shear. This is true also 
for time dependent shears, see |2T| 120] and references therein. For each 
realization, c*(S, u) = c* Q + 0(5 2 ) as 5 < 1; c*(5,u) = c* + 0(5) as 5 > 1. 
For each realization, as 5 <C 1, we have: 

Proposition 2.1 Let \ = x(y) solve A y x = —b, y G Q, with zero Neumann 
boundary condition, where b G C(Q), has zero mean over Q. Then for 5 
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sufficiently small, the minimal speed has the expansion 



c*(6) = Co + ^ / I V X | 2 dy + 0(5 3 ). (2.5) 

I I J fa 

Up to 0(5 2 ), the formula (|2.5|) is independent of the nonlinearity, see [221 
for / being bistable or combustion nonlinearity, also [21] for more general 
nonlinearities and time periodic shears. We shall give a proof of Proposition 
2.1 using variational formulas, and later generalize it to the random case. A 
helpful fact is that the infimum in (jl.2)) can be restricted to a bounded set 
independent of b and S, as stated in the following lemma: 

Lemma 2.1 Let b G C(0) have zero mean over tt, and let A = y/f'(0). 
Then 

inf^= inf (2.6) 

A>0 A 0<A<A A 

Proof of lemma: For each c > 0, we let p c (A) = /i(A) — Ac — A 2 . So, if cf> > 
is the eigenfunction defined by (jl.Hj) . then p c (A) is the principal eigenvalue 
defined by the equation 

A y 4> + [Xb(y) - Ac + /'(O)]0 = p c (A)0, t/6fi. (2.7) 

One can readily verify that d\ p c (A)|a=o = — c < 0. The variational formula 
(ll.2j) can be expressed as 

c* = inf {c | 3A > 0, Ac = /i(A)} 

= inf {c | 3A > 0,p c (A) = -A 2 } . (2.8) 

Consider the points where p c (A) = —A 2 . By Proposition 2.1 of [§], the 
continuous curve A i— > p c (A) is convex in A, for each c > 0. Also, p c (0) = 
/'(0) > 0. Therefore, for a given c > 0, there can be at most two values of 
A > such that p c (A) = -A 2 . The line p*(A) = -2y/f (0)A + /'(0) satisfies 
P*(A) > —A 2 , with equality holding only at one point: A = a/ /'(0). Since 
p*(0) = p c (0) and p c (A) is convex and p* is a line, p c (A) = —A 2 for some 
A > only if p c (Ai) = —A 2 for some Ai G (0, A ]. This point is illustrated 
in Figure ^ The solid curve represents the parabola —A 2 . If p c (A) intersects 
—A 2 , then one of the intersection points must be to the left of Ao = a/ f'(0). 
Therefore, from (J2.8)) . 

{c | 3A > 0,p c (A) = -A 2 } = {c | 3A G (0, A ],p c (A) = -A 2 } . (2.9) 
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Figure 1: Intersecting curves p c and —A 2 . 



So, we conclude that 



c*(5) = M^= inf (2.10) 

0<A A 0<A<A(, A 



Proof of Proposition 2.1: 

To estimate c*(S), we bound the principal eigenvalue //(A) using two dif- 
ferent representations of \i. First, since L is a self- adjoint operator, we have: 



where the supremum is taken over all V e # 2 (^) such that §£ = on dQ. 
The other representation is: 

H = inf sup ^ = inf sup ( ^f- + A56 + A 2 + /'(0)^| , (2.12) 

<P yen W <Pyen\W J 

where the infimum can be taken over all ip G C 1 (f2) such that Aip e C(f2), 
■0 > 0, and |^ = on <9f2. This representation follows from the fact that 
the eigenfunction <fi > lies in the kernel of the self-adjoint operator {L\ — 
/x(A)J) = (L\ — /i(X)I)*. So, if we have the strict inequality 

L x ip - = m < (2.13) 
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then the Fredholm alternative implies that {(f), m)^ = 0, a contradiction 
since <fi > 0, m < 0. Hence, 

sup ^> MA). (2-14) 

yen W 

Since L\(f> = fi(X)4>, the formula ()2.12j) follows. Note that we do not require 
the test functions ip to be C 2 (Q), only Aif> e C(Q). This is important since we 
do not want to require the shear b(y) to be any more regular that b G C(Cl). 

Let us derive upper and lower bounds for MA) by choosing test functions 
if) as: 

rj} = 1 + \5 X + \ 2 5 2 h, (2.15) 
where x — x(y) an d h = h(y) solve 

A X = -b 

Ah = -bx + k, (2.16) 
with zero Neumann boundary conditions at dfl, and k a constant equal to: 

k = ^ f ' b Xdy = j^ [ |Vx| 2 dy. (2.17) 



We normalize x an d h so that: 



infx(») = and inf /i(x) = 0. (2.18) 



Then 

L A ^ = A 2 <5 2 A; + A 3 5 3 6/i + (A 2 + /' (0))i/j, 

and 



, - X 2 S 2 k^- + X 3 S 3 ^f + A 2 + f (0). (2.19) 

2 W )W 

2 _ _ \2r2 



Using the definition of if), we see that if) = ip — \5x4> ~ X 5 hip, and 

L^ 2 dy L^ 2 dy L^dy 

Now from (j2.11|) and ()2.19|) we have the lower bound 

MA) > X 2 + f'{0) + X 2 5 2 k + R l , (2.20) 
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with 

i?i = - J"^ (k j XiJ + k\5 J hif) - J bhif^J (2.21) 

By choice of x > and /i > 0, we have j n ip 2 > |^| for all 5 > and A > 0. 
Hence, i?i = 0(<5 3 ) for A bounded. Returning to the variational formula 
(I2.10JI . we now have a lower bound on c*(S): 

c*(6)= inf 4^ > inf fA + ^ + A^ + 4 1 

0<A<A A 0<A<A \ A A 

> inf ( A + ^ + A5 2 ^ + 0(<5 3 ) 
~ a>o \ A J 

= 2y/f'(0)(l + 6 2 k) + O(5 3 ) 

= c + ^p + O(5 3 ). (2.22) 

To obtain an upper bound on c*(S), we use f!2.12|) and calculate: 
L ^ A ^+X5b + X 2 + f'(0) 



if) ip 

X 2 5 2 k + X 3 5 3 bh 



1 + \5 X + X 2 S 2 h 
Since x > and h > 0, we see from (j2~T2|) and (j2~28|) that 



+ A 2 + f(0). (2.23) 



At (A) < sup ^ < A 2 + f(0) + A 2 <5 2 £; + i? 2 , (2.24) 

yen V 

with 

i? 2 = A^Hft/ilU. (2.25) 
The variational formula (j2.10|) implies: 



c*(5) = inf ^-^ < inf A + + A<Tife + R 2 

0<A<A A 0<A<A \ A 

= c + ^p + O(5 3 ), (2.26) 

completing the proof. 

When the shear b(y, u) is a random process, the corresponding minimal 
speed c*(S) = c*(S,u) is a random variable for each 5, and we consider how 
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the expectation £"[c*(5)] scales with the parameter 5 by finding an exponent 
p such that J5[c*(5)] = c*(0) + 0(5 P ). Each realization of the process b(y,u) 
restricted to the domain Q does not necessarily have zero integral over Q. 
Nevertheless, each realization can be written in the form 

b{y,u) = b{u) + b 1 (y,Lo), (2.27) 

where b(u) = (b(y,u)) is the mean of b over Q, and b\(y,u) is the variation 
about the mean value. For a fixed realization, the minimal speed c*(S) will be 
affected by both the scaling of the mean b(u) and the scaling of the varation 
bi(y, uS). That is, 

c*(6,u) = (%-6b(u) + M(6,w), (2.28) 

where the remainder M(8, u>) is the enhancement due to the variation b\(y, u), 
different for each realization. Taking the expectation of both sides of (|2.28|) . 
we have 

E[c*(5)} =c* - 6E[b(u)} + E[M(6, u)] (2.29) 

For each sample, M{5,u) is 0(5 2 ) for 5 small. Though E[M{$)\ might ex- 
hibit different scaling than quadratic, we show that the quadratic scaling law 
remains for enhancement of averaged front speeds under suitable moment 
conditions of the shear. 

Theorem 2.1 Let b(y,u) be a stationary random process in R 71 ^ 1 (n > 2) 
so that sample paths are almost surely continuous; and that 

E[\\b\\l] < +oo. (2.30) 

Then for 5 small, the expectation E[c*(8)] has the expansion 

E{c*(5)} = co - 5E[(b)} + ^-J n E{\VX\ 2 } dy + 0(5 3 ), (2.31) 

where b(y, u>) = (b)(u) + bi(y, uj); and x = x{Vi solves A y x = — 6i, y & ; 
subject to zero Neumann boundary condition. 

Proof: As the contribution of {b) to c* is just an additive constant, it suffices 
to consider shear flow b\ and show that it gives the averaged speed 

E[c*(S)} = c + ^ f E[\V X \ 2 ] dy + 0(<5 3 ). (2.32) 

I | «/ Sa 
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We adapt the proof of Proposition 2.1, noting however that in the stochastic 
case the remainders R\ and R 2 defined by (|2.21j) and (J2.25|) are random and 
not bounded uniformly for all realizations. Instead, we will show that for A 
in a bounded interval, 

E[\Rx\) < 0(5 3 ) and E[\R 2 \] < 0(5 3 ). 

To this end, we estimate x an d h, with C denoting a generic positive constant 
depending only on the domain Q and its dimension. Let x an d h solve (|2.16|) 
with (x) = (h) = 0. Applying W 2,p estimates (e.g. Thm. 19.1 of 9J), we 
have 

<CM 1/p ||&i||co (2.33) 



and 



w 2 'P{n) < C\\bix + k\\ LP (u) 

< CH^IIooIIxIIl^+CA;^! 1 ^ 

< CWhWl (2.34) 



since 



fc = (|Vx| 2 ><C||&X (2.35) 

Given a G (0,1), we can choose p > 1 sufficiently large, depending on n, 
such that W 2 ' P (Q) embeds continuously into C 1,a (Cl). It follows that there is 
a constant C > independent of b such that 

Hxllci(fi) < C||6i||oo and ||/i|| C i(n) < C||^i||L ( 2 - 36 ) 

If instead we normalize x an d h by ()2.18|) . then the bounds ()2.36|) still hold, 
with different constants. Note that adding a constant to x an d h does not 
alter the quantity J Q |Vx| 2 that appears in the asymptotic expansion. 
Now by (|2.36|) . the integrals in R 1 are easily bounded as: 



xi> = x + ^x 2 + ^s 2 xh 

Jn 

< C'(||6 1 || 00 + A5||6 1 || 2 <) + A 2 5 2 ||6 1 ||L) 



Similarly, 



hip = J h + X5 X h + X 2 5 2 h 2 
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and 



bihip dy 



b x h + X5 j bih X + X 2 S 2 J b x h 2 
< C(||6 1 || 3 00 + A5||6 1 || 4 00 + A 2 5 2 ||6 1 || 5 00 ). 



Since \ an d h are nonnegative, f n ip 2 dy > C f n ip = C\Q\ > 0, for any 
realization. So for A in a bounded interval and 5 small, we bound (|2.21|) by: 

l^il < C<5 3 A 3 (l + ll&iHD, 

so that £[|#i|] < 0(5 3 ). 

To bound £/[|i? 2 |], we use the normalization (|2.18|) and the above esti- 
mates: 



A 3 5 3 \\bxh\U < CA 3 5 3 ||&i 



3 



|3 

I oo' 



(2.37) 



Hence E[R<2\ < 0(5 3 ) for A in finite interval. Now we return to ()2.22|) to 
conclude 

E[c*(5)] > E[2^f'(Q)(l + 6 2 k)} + 0(5 3 ) 



Co 



c 5 2 E[k] 



+ 0(5 3 



since E[k 2 } < CE[||6||^] < oo, 

The opposite inequality follows from ()2.26|) since E[R 2 ] = 0(5 3 ) for A G 
(0, Aq). Thus formula (|2.32|) holds. For general b, not necessarily mean zero, 



E[c*(6)} 



c + 6E[(b}} + 
c + 6E[(b}} + 



col 
2 

2\n\ 



E[k] + 0(5 3 ) 

E[\W X \ 2 ] dy + 0(5 3 ) 



The proof is finished. 

In our numerical computation on front speeds in two dimensional chan- 
nels, we use the Ornstein-Uhlenbeck (O-U) process for shear b, so E[(b)] = 0. 
Let us show below that the O-U process, denoted by X(y,u), satisfies the 
conditions in Theorem 2.1, and so S[c*(5)] scales quadratically with 5 for 5 
small. 
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Corollary 2.1 (Explicit Average Speed Formula) Consider the O-U 
process b(y, uS) as solution of the Ito equation: 

dX(y) = -aX(y)dy+rdW(y) } ye[Q,L], (2.38) 

where W{y,uS) is the standard Wiener process, X(0,u) = X (uj) is a Gaus- 
sian random variable with mean zero, and variance p = r 2 /(2a). Then 
X(y,uj) satisfies the moment conditions in Theorem 2.1. The averaged KPP 
front speed in the channel R x [0, L] is given by 

E[c*(5)]=c + ^enh + O(5 3 ), 8 < 1, (2.39) 

where: 

r 2 ( _ aL ( 4 1 \ L 4 5 4 



enh = — e 



2a V \L 2 a A 3a 2 J 3a L 2 a 4 3a 2 La 3 

Proof: The O-U process is stationary and Markov. Its sample paths are 
almost surely Holder continuous though nowhere differentible. The process 
can be written as 

b(y, to) = e~ ay b(0, u) + r f e~ a{y ~ s) dW s (u), (2.40) 



The covariance function of this process is pe a ' y S L Letting g(y,u>) denote 
the process 

g(y) = e a yb(y, u) = g(0, u) + r f e as dW s (u), (2.41) 

we see that g(y,u) is a martingale [T3*] . By Doob's martingale moment 
inequality [TH], for any p G (1, +oo), 

E[ sup \g(y)\*\ < f-P-X E[\g(L)n (2.42) 

0<y<L \p-LJ 

Since the process b(y, lu) is Gaussian, (|2.41jl and ()2.42j) imply that 

E[\\b\t]<CE[\b(L)f}<+oo. (2.43) 
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Formula (|2.39|) now applies to the average speed. Notice that 

bi(s)bx(y) ds dy, 

and 

ElhisMy)} ds dy (2.44) 
Let us calculate E[bi(s)bi(y)} in terms of E[b(s)b(y)]. Define: 
g(y) = (/(•, y)> or g(s) = (f(s, •)) 

so that 

EMyMs)} = E[b{s)b{y)\ - E[b{s)b\ - E[b(y)b] + E[b% 

1 ' L 



E[b(s)b} = - I E[b(y)b(s)]dy = g(s), 



o 



E[b 2 ] = j- 2 / / E[b(s)b(y)]dyds=(g). 



o Jo 



Thus 

EMvMs)] = f(s, y) + (g) - g(y) - g(s). 

Now, we have 



x rx 



E[(Xx(x)Y} = / / ElhisMy^dsdy 
Jo Jo 



x rx 




f(s, y) + (9) - g(y) - g(s) ds dy 

rx rx 

x 2 (g) - 2x 2 (g) x + / / f(s,y)dsdy, 



'0 jo 



JO 



where (g) x denotes the average of g over the interval [0,x], for < x < L. 
Consequently, we have 



1 " L 



E[(\Xx\ 2 )} = j I [x 2 (g)-2x 2 (g) x + 







nx \ 
f(s,y) ds dy) dx 
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Using the 0-U covariance function, we proceed as: 

»w = \[ /(•.»)*= si f e " l "~' ,dy 

r 2 f \ — g _as x — e - a (L-s) 

+ 



and 



2a \ La La 



r 2 1 fV-e^ 1 -e^ - ") , 

<7* = 7T- / (— ? + F ) ds 

la x l n La La 



r 2 / 2 1 1 

2a Ua xLa 21 j + xLa 2 1 j 



Letting x = L, we have 



Similarly, 




'0 </o 

Combining the above, we have 



r 2 f It 2 
f(s,y) ds dy = -( — + -(e~ ax - 1) ) . 
2a \ a cr 



r 2 2 [ L 2x 2 x , _ _ x 



, , , - +-^-( e - a;!; -l) + — (e- ai -e- a ( L -^)^ 
2a L /n La La z La z 



o 

r 2 1 f L 2x 2 



+—- / — + -( e — 
2a L J a a A 

r 2 ( _„ T ( 4 1\ L 



e aL 



2a V V £ 2 a 4 3a 2 / 3a L 2 a 4 3a 2 La 3 / ' 

(2.45) 



In view of (I2.39|) . the proof is complete. 
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Theorem 2.2 (Linear Growth) If the stationary shear process b(y,u) has 
almost surely continuous sample paths and satisfies ElWbW^} < oo, then the 
amplified shear field 5b(y,u>) generates the average front speed: 

E[\c*(6,w)\]=0{6), <5>1. 

Moreover, lim^oo .E[|c*((!>, u))\]/8 exists. 

Proof: By Theorem 5.1 of [4 , l c *^'")l — > d*(u) as 5 — ► oo, where d* is finite 
for each to. Now recall the upper bound \c*(S, < |c | + 5||6||oo- Hence 
for 5 > | c U < 1 + IHL = Y, and E(Y) < oo. The dominated 

convergence theorem implies that: 

E[^^]^E[d*(u)}<E(Y). 
o 

The proof is finished. Clearly, the O-U process satisfies the required condition 
for linear average speed growth. 



3 Computation by Variational Principle 

3.1 Numerical Methods 

Let n = 2. For a given A > 0, we compute the principal eigenvalue //(A) with 
corresponding eigenfucntion (f> = <p(y) > 0, y G [0, L], by solving: 

<P yy + [A 2 + A b(y) + f(O)]0 = /i(A)0, y e (0, L), 

^ = 0, y = 0, L, (3.1) 

using a standard second order finite-difference method. Here we suppress 
the random parameter uj, as computation is done realization by realization. 
Denote the uniform partition of the domain by points and the nu- 

merical solution by (f> = {<^i}._-., where h — L/(m— 1), y% = (i — l)h, and 
(pi m <p{yi)- The discretized system is 

1 2 - 1 - 

— + (A 2 + \h + f (0) - — )0i + — m = /i(A)0i i = 2, . . . , m - 1, 

with second order approximation of the Neumann boundary conditions. This 
reduces to finding the principal eigenvalue of a symmetric tridiagonal matrix, 
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Figure 2: Two curves A i— > 

easily accomplished with double precision LAPACK routines [2|. Then we 
compute points on the curve H(\) = and minimize over A using a 

Newton's method with line search. Our approximation decreases with each 
iteration and converges quadratically in the region near the infimum. Two 
illustrative curves H(X) are shown in Figure 01 for two different realizations 
of the shear. 

We generate realizations of the shear process b(y, u) by numerically eval- 
uating the stochastic ODE ()2.38|) with the Milstein scheme (see [IH])- Al- 
though this scheme is first order, we use a discrete spacing h < h 2 , where h 
is the discrete grid spacing for the eigenvalue problem, so that the method is 
still second order accurate in the parameter h. Figure El shows a sample path, 
and Figure E] compares exact and numerical covariance function constructed 
from 5000 samples. 

To approximate the expectation £"[<;* (5)] we generate N independent re- 
alizations (indexed by i = 1,...,N) of the shear and compute the corre- 
sponding minimal speeds {c*} for each 5. Then we compute the average 
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Figure 3: One sample path of the Ornstein-Uhlenbeck process b(y,u). 




Figure 4: Numerical and exact covariance functions of the Ornstein- 
Uhlenbeck process b(y,u). 



E[c*(5)} « E{8) = c* + i £ M t (5), (3.2) 

t=i 

where Mj(5) = c*(<5) — Cq + That is, we subtract the linear part due to 
the mean of the shear being nonzero, as in (|2.28jl . 

Once we have the averages E(5) for each 5, we compute the exponents p 
using the least squares method to fit a line to a log-log plot of speed versus 
amplitude. That is, the exponent p is the slope of the best-fit line through 
the data points (log(5), \og(E[c*(5)} — c*(0))) for each shear amplidute 5. 
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Figure 5: Average enhancement of minimal speed in small amplitude shears. 
The solid curves are the numerically computed values; the dashed curves are 
given by formula (|2.39Jl . 

3.2 Numerical Results 

3.2.1 Scaling with shear amplitude 5 

In Figure 03 and Figure |H1 we show the results of a simulation using iV = 
100, 000 realizations of a shear in small and large root mean square (rms) 
amplitudes, respectively. As shown later in Figure ITT1 we find that the choice 
of N = 100, 000 samples was more than enough to obtain good convergence 
of the speed distribution functions. The covariance function of the process is 
E[b(y)b(s)] = 2e~ 4 ' y ~ s '. In each plot, we show multiple curves, corresponding 
to various domain sizes. In Figure El corresponding to small 5, the solid 
curves are the numerically computed values; the dashed curves are given by 
formula (|2.39jl . We find that the enhancement of the minimal speed scales 
quadratically for small amplitudes and linearly for large amplitudes. The 
computed exponents are shown in Tabled 

As an application of the numerical results for large amplitudes, we com- 



17 



Table 1: Computed scaling exponents p for i£[c*(5)] = Cq + 0(5 P ). 
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Figure 6: Average enhancement of minimal speed, large amplitude shears. 

pare the distribution of c*{8) with the distributions of the random variables 

gi{uj) = 2k/'(0) + 5||6||oo and ^(^) = 2yK + ^|jVx||oo, where k is the dif- 
fusion constant (equal to 1 in equation Theorem 2 of [TT] shows the 
upper bound: 

c*(6,u) < m.ia(g 1 (u),g 2 (uj)) (3.1) 

and #1 < g>2 provided that ^jtj^ is sufficiently small, depending on the 

realization. In Figure [3 we compare the distributions for c*(50), g±, and g 2 
for k = 0.01. The asymmetry is seen in all three curves. 
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Figure 7: Distributions of speed enhancement, lower and upper bounds at 
5 = 50. 

3.2.2 Dependence on Covariance 

Next, we consider the effect of the covariance on the enhancement of the 
minimal speed. The covariance E[b(y)b(s)] is a function of \t\ — \s — y\, so 
we will write V(t) = E[b(y)b(s)]. By choosing r = v^2« 3//4 , we constructed 
O-U processes with covariances given by 

V(t) = V^e~ altl . (3.2) 

By this choice of r, the L 2 norm of V(t) remains constant as a changes, so 
that the total energy in the power spectrum of the signal remains constant. 
Since |- = yfa, we see from equation ()2.39j) that for fixed L, 

lim E[(\ Xx \ 2 )]= lim . E{(\ Xx \ 2 )]=0 (3.3) 

and that i?[(|x x | 2 )] achieves a maximum for some finite value of a G (0, oo). 
This suggests that there is some optimal a, depending on the domain size L, 
such that the enhancement of £J[c*(5)] is maximized. 

Fixing the grid spacing dx = 0.002, we computed the expected value 
i£[c*(5)] for a range of a and for L = 1.0,2.0,3.0,4.0. Note that for each a, 
we must choose the initial points b to have variance E[bo\ = \fa so that the 
process remains stationary for each a. Varying the covariance does not effect 
the order of the scaling in 5. That is, in each case the enhancement scales 
like 0(5 2 ) for small S and 0(5) for large 5, as in the preceding simulation. 

Figure |H1 shows the enhancement £'[c*(5)] for a fixed 5 = 1.0 and a range 
of a. Figure El shows the results of the same computation for 5 = 15.0, 
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Figure 8: Effect of covariance on minimal speed enhancement at 5 = 1.0. 

corresponding to the large S regime. In this case, formula (|2.<fflj) is no longer 
valid. Nevertheless, we see the same effect as in the small amplitude regime: 
the existence of an optimal covariance parameter a. 

This effect can be interpreted in terms of V(t) and its Fourier transform 
or power spectrum: 



V(w) = 




As a — > 0, V(w) concentrates at the origin, and so the energy of the shear 
process is concentrated more in the large scale spatial modes. The domain 
Q, to which the process is restricted, is bounded, and variations over a length 
scale that is much greater than the diameter of Q have little effect on the 
average enhancement of the front. As a result, E[c*] decreases as a — > + . 
In the other limit a — > oo, V spreads out so that the energy over any finite 
band of frequencies goes to zero, causing E[c*} to decrease as well. Note 
that V — > in L l as a —>■ oo, so even though V spreads out more uniformly 
as a — > 0, the family of processes does not converge to white noise, whose 
covariance function equals the Dirac delta function. 
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Figure 9: Effect of covariance on minimal speed enhancement at 5 = 15.0. 
3.2.3 Speed Distribution 

For a fixed 6=1 and 5 = 14 (corresponding to small and large amplitudes), 
we computed the distributions of the numerically computed values M{$). 
The distributions are shown in Figure To compute these distributions, 
we partition the range of values into Q disjoint intervals: {[xj, Xj+i)}- =r 
Then, we let 

Pdf{x) = ^^^- ) if xe[ Xj ,x j+1 ) (3.5) 

where Xj( x ) is the characteristic function of the interval [xj, Xj+i). The 
distributions in Figure HO were computed with N = 100, 000 samples, and 
Q = 300. 

The values M(5) are the enhancement of the minimal speeds due to the 
variation of the shear, after the effect of the mean field has been subtracted 
off. Since a mean zero shear always enhances the minimal speeds, we should 
expect M{$) > for all 5, for all realizations. As already noted, Figure 
ITT1 shows good convergence of the distributions when using N = 100, 000 
samples. 
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Figure 10: Computed probability distribution functions (pdf's) of enhance- 
ment M(5), S = 1.0 (bottom), 5 = 14.0 (top). 



3.2.4 Comparison with Direct Simulation 

When comparing the results of the previous sections with results from direct 
simulation, we find that computing the minimal speeds using the variational 
formula offers significant advantages. As in |21j . we first computed the quan- 
tities £'[c*((5)] via direct simulation of the original equation (jl.ljl on a trun- 
cated domain using an explicit second-order upwind finite difference scheme. 
We chose the diffusion constant in the direct simulations to be k = 0.025 
and the grid spacing to be dx = 0.05, dt = 0.004. For each realization of 
the shear, we evolved the solution of (jl.ljl for a sufficiently long time until 
the front moves at a more or less constant speed (see j2l] for details of the 
method). In this way, we approximated the minimal speed for each realiza- 
tion; then, we repeated the process for a large ensemble of shears to compute 
the expectations E[c*(S)]. 

We considered three nonolinearities in our direct simulations: 

• KPP nonlinearity: f(u) = u(l — u). 

• Combustion nonlinearity: f(u) = for u G [0,8], and f(u) > for 
u e (6, 1], for some 6 e (0, 1). Also, f (1) < 0. 
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Figure 11: Convergence of speed enhancement distribution at 5 = 14.0. 

• Bistable nonlinearity: f(u) = u(l — u)(u — fx) for some G (0,1/2). 

Although there is no known variational formula as simple as formula ()1.2|) 
for the bistable and combustion nonlinearities, the expansion ()2.5|) holds for 
each of the nonlinearites (see also Theorem 4.2 of ^2j). Therefore, we should 
expect that for small 5, the computed values 

^«^/|Vx(!/,-)| 2 A/ + OOT (3.6) 

have approximately the same distribution, independent of the nonlinearity. 
For each nonlinearity, we computed more than 700 realizations of the shear 
and evolved the solution, as described in |21j . A relatively smaller number 
of samples was due to the long time required to compute each sample (not 
a problem with the variational formula in the KPP case). After subtracting 
off the linear part of the enhancement (due to mean field), we computed the 
enhancement exponent as described in Section 3. The results of the direct 
simulation concur that, for small amplitudes, _E[c*(<5)] scales like 0(6 2 ). Fig- 
ures[T21[ini and ITU show the distributions of the computed values for the 

c 

KPP, combustion, and bistable nonlinearities, respectively, for L — 1. We 
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see that the distributions are roughly the same, independent of the nonlin- 
earity, as should be expected. While computationally expensive, the direct 
simulation method reveals the universal 0(S 2 ) scaling of the front speeds for 
small 5 (albeit on smaller ensembles than in the variational computations), 
for each of the nonlinearities considered. 

Because the shears are random and may greatly distort the wave front, 
one major challenge in accurately approximating the (minimal) speeds is 
tracking the widely varying front region over a long time. In the region 
of the front, gradients are relatively large, so accurately tracking the front 
requires either a very fine uniform grid spanning a large domain or some 
kind of adaptive-mesh scheme. In either case, accurate direct simulation 
is prohibitively expensive compared to the simple variational method. In 
contrast to the direct simulation method, the variational formula allowed us 
to compute a much larger number of samples in a fraction of the time and 
avoid the inaccuracies resulting from truncation of the channel doman. 

Finally, we note that if Q is unbounded, then the quadratic asymptotic 
behavior of c*(S) cannot hold in general. For example, in it was shown 
that if the channel R x [0, L] is replaced by R 2 then the front velocities 
obtained through a comparable variational principle diverge due to the al- 
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Figure 13: Distribution of speed enhancement via direct simulation, combus- 
tion nonlinear ity, 5 = 0.5. 

most sure growth of the running maximum of the process b(y, u). This effect 
can be clearly seen in our results, as i?[(|Vx| 2 )] given by (|2.39|) diverges 
as L — > +oo. In our application of Theorem 2.1 to a stationary Gaussian 
process, the boundedness of the domain Q is crucial in order to achieve the 
necessary bound on ||6||oc (through the martingale moment inequality), and 
in turn the bound on c*. 



4 Conclusions 

Sufficient moment conditions are obtained to ensure the quadratic (linear) 
KPP average front speed enhancement through small (large) rms random 
shear flows in channel domains. The conditions are realized by the Ornstein- 
Uhlenbeck process for which an explicit average front speed formula is de- 
rived. The variational principle based computation is carried out for the 
speed ensemble. Numerically computed speed enhancement is in agreement 
with theoretical analysis, and provides data for studying speed distributions 
and dependence on shear covariance. Comparison with direct simulations of 
random fronts in case of non-KPP nonlinearities showed the same enhance- 
ment scaling laws. It is interesting to investigate front speeds through time 
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Figure 14: Distribution of speed enhancement via direct simulation, bistable 
nonlinearity, 5 = 0.5. 

dependent random shear flows or non-shear random flows in channel domains 
for future studies. 
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